Introduction

In this report I will be finalizing the LME model for Padilla zebrafish LMR assay for DNT Loperamide exposure data. Light data model was never finalized. Data was never fit with centered time data. In Light there was also signs in the residuals that the inclusion of a cubic polynomial time may benefit the model. These two operations will be performed. Then Dark data will be fit, CIs for Light and Dark sample level parameters will be graphed, and significance tests evaluated.

Load and Visualize Data

Loading and transformation of data will not be shown. Time data will be centered and movement data will be square root transformed. Light and Dark data will be plotted separately.

Plot individual trends, Light and then Dark.

Fit Models to Light and Dark Data

Light data was never evaluated with centered time. Evaluate the centered time data Light models. Fit cubic model to see if it is preferred by AIC.

ctrl <- lmeControl(opt = "optim", niterEM = 200, msTol = 1e-20, msMaxIter = 1000, apVar = FALSE)

lmm1.L <- lme(sqrty ~ conc*poly(t,degree=2,raw=TRUE),
            data = Light.sqrt,
            random = reStruct(~1|fishID,REML=FALSE,pdClass="pdSymm"),
            correlation = corAR1(,form=~t|fishID),
            method = "ML",
            control = ctrl)
lmm2.L <- lme(sqrty ~ conc*poly(t,degree=2,raw=TRUE),
            data = Light.sqrt,
            random = reStruct(~poly(t,degree=1,raw=TRUE)|fishID,REML=FALSE,pdClass="pdSymm"),
            correlation = corAR1(,form=~t|fishID),
            method = "ML",
            control = ctrl)
lmm3.L <- lme(sqrty ~ conc*poly(t,degree=2,raw=TRUE),
            data = Light.sqrt,
            random = reStruct(~poly(t,degree=2,raw=TRUE)|fishID,REML=FALSE,pdClass="pdSymm"),
            correlation = corAR1(,form=~t|fishID),
            method = "ML",
            control = ctrl)
lmm4.L <- lme(sqrty ~ conc*poly(t,degree=3,raw=TRUE),
              data = Light.sqrt,
              random = reStruct(~poly(t,degree=3,raw=TRUE)|fishID,REML=FALSE,pdClass="pdSymm"),
              correlation = corAR1(,form=~t|fishID),
              method = "ML",
              control = ctrl)

zz <- AIC(lmm1.L,lmm2.L,lmm3.L,lmm4.L)
normAIC(zz)
##        df       AIC
## lmm4.L 52  0.000000
## lmm2.L 35  4.392619
## lmm3.L 38  7.295060
## lmm1.L 33 22.582679
## Random effects:
##  Formula: ~poly(t, degree = 1, raw = TRUE) | fishID
##  Structure: General positive-definite
##                                 StdDev     Corr  
## (Intercept)                     0.55419666 (Intr)
## poly(t, degree = 1, raw = TRUE) 0.02385692 1     
## Residual                        1.00000000
## Random effects:
##  Formula: ~poly(t, degree = 2, raw = TRUE) | fishID
##  Structure: General positive-definite
##                                  StdDev      Corr          
## (Intercept)                      0.624589806 (Intr) pd=2r=T
## poly(t, degree = 2, raw = TRUE)1 0.022983985  1            
## poly(t, degree = 2, raw = TRUE)2 0.001478449 -1     -1     
## Residual                         1.000000000
## Random effects:
##  Formula: ~poly(t, degree = 3, raw = TRUE) | fishID
##  Structure: General positive-definite
##                                  StdDev      Corr                  
## (Intercept)                      0.657772171 (Intr) p(,d=3,r=TRUE)1
## poly(t, degree = 3, raw = TRUE)1 0.057159545 -0.244                
## poly(t, degree = 3, raw = TRUE)2 0.002295310 -0.696 -0.055         
## poly(t, degree = 3, raw = TRUE)3 0.000800015  0.619 -0.891         
## Residual                         1.000000000                       
##                                                 
## (Intercept)                      p(,d=3,r=TRUE)2
## poly(t, degree = 3, raw = TRUE)1                
## poly(t, degree = 3, raw = TRUE)2                
## poly(t, degree = 3, raw = TRUE)3 -0.128         
## Residual

Amazing how much correlation of random effects changes with the introduction of a cubic term. Visualize this behavior.

We need to evaluate the residuals of the models and the accuracy of parameter estimates to see if using the cubic model in Light would actually be beneficial for the purposes of evaluating chemical effect.

Evaluate Models

First, evaluate parameter estimation accuracy.

Confidence around slope parameter estimates is much smaller for lmm2.L. The introduction of 10 more parameters, a cubic parameter for every concentration group, may explain the loss in accuracy. For our purposes lmm2.L model may be more ideal. Evaluate the change in parameter accuracy b/t lmm2.L and lmm3.L. Then visualize the fits and residuals.

Little to no difference except at intercept term. Visualize the fits.

Compare Fitted models

Visually evaluate the fitted models for lmm2.L, lmm3.L, and lmm4.L. Start with lmm2 vs. lmm3. nlme plotting functions limit us to 2 models at a time.

Fits are essentially identical.

Cubic term adds slight cubic appearance to the fits but the fits are super similar. For some subjects, not observed in this random sample, the cubic term may be more beneficial.

Compare the population-averaged models.

Cubic term doesn’tadd curvature at the beginning of the Light phase for all concentration groups. I was hoping it would do this. Cubic function appears to model some concentrations better than others.

Plot Model Resiudals

Evaluate normality of residuals, residuals by time and residuals by concentration.

Residual values look okay for sample level (level = 0), but you can see at the tail end what is likely the effect of zero values in the Light data. Comparing lmm2.L and lmm4.L there is surprisingly little to no change in the sample-level residuals. At the level of the subject it seems like fits are better. I am going to stick with lmm2.L.

Before I finish this section I want to mention one problem of these models, both Light and Dark, the heterogeneity of time periods. In the nlme book I purchased there are variance functions described that can be applied to deal with this heterogeneity. One in particular can model variance as a power function of time. I will show the residuals of lmm2.L against time to demonstrate this power-like behavior. The addition of this variance function would increase the number of parameters estimated, so the benefit of doing so is questionable. In a further iteration of this model I imagine this being something we attempt to do.

Refit lmm2.L with REML. This will be our final Light phase model.

lmm2.L.REML <- lme(sqrty ~ conc*poly(t,degree=2,raw=TRUE),
            data = Light.sqrt,
            random = reStruct(~poly(t,degree=1,raw=TRUE)|fishID,REML=TRUE,pdClass="pdSymm"),
            correlation = corAR1(,form=~t|fishID),
            method = "REML",
            control = ctrl)

Visualize the population fits by concentration group.

Dark Model

Fit Dark model with REML. Dark model identified in Week 04-11-2022 report will be used. Evaluation of that model was performed in that report. Note that “lmm3” full quadratic model will be used. I have yet to perform Aikake weight calculations.

ctrl <- lmeControl(opt = "optim", niterEM = 200, msTol = 1e-20, msMaxIter = 1000, apVar = FALSE)

lmm3.D.REML <- lme(sqrty ~ conc*poly(t,degree=2,raw=TRUE),
            data = Dark.sqrt,
            random = reStruct(~poly(t,degree=2,raw=TRUE)|fishID,REML=FALSE,pdClass="pdSymm"),
            correlation = corAR1(,form=~t|fishID),
            method = "REML",
            control = ctrl)

Visualize the population fits by concentration group.

Overall Dark model appears to fit the Dark phase data much better than Light model fits Light.

Plot Confidence Intervals in Conc-Resp Mode and Evaluate t-Tests

Light parameters.

Dark parameters.

Evaluate t-tests for Light.

##                       Value   Std.Error   DF     t-value p-value
## conc0.004 B_0  1.027109e-01 0.207235857  230  0.49562299   0.620
## conc0.012 B_0  4.369704e-01 0.207235857  230  2.10856552   0.036
## conc0.04 B_0  -1.881950e-01 0.207235857  230 -0.90812009   0.364
## conc0.12 B_0  -9.487428e-02 0.207235857  230 -0.45780821   0.647
## conc0.4 B_0   -1.504684e-02 0.207235857  230 -0.07260732   0.942
## conc1.2 B_0    1.830925e-02 0.207235857  230  0.08834980   0.929
## conc4 B_0     -3.813707e-01 0.227015507  230 -1.67993242   0.094
## conc12 B_0    -5.484878e-01 0.227015507  230 -2.41608056   0.016
## conc40 B_0    -1.081224e+00 0.219690489  230 -4.92157956   0.000
## conc0.004 B_1 -5.748389e-03 0.017002290 4540 -0.33809496   0.735
## conc0.012 B_1 -1.400175e-02 0.017002290 4540 -0.82352152   0.410
## conc0.04 B_1   1.284161e-03 0.017002290 4540  0.07552869   0.939
## conc0.12 B_1  -2.926721e-02 0.017002290 4540 -1.72136862   0.085
## conc0.4 B_1   -5.582097e-04 0.017002290 4540 -0.03283144   0.973
## conc1.2 B_1   -2.258041e-02 0.017002290 4540 -1.32808069   0.184
## conc4 B_1      2.901809e-04 0.018625076 4540  0.01558012   0.987
## conc12 B_1    -1.792078e-02 0.018625076 4540 -0.96218548   0.336
## conc40 B_1    -4.308857e-02 0.018024108 4540 -2.39060738   0.016
## conc0.004 B_2 -2.629120e-03 0.002641350 4540 -0.99536979   0.319
## conc0.012 B_2 -2.719347e-03 0.002641350 4540 -1.02952925   0.303
## conc0.04 B_2   4.549159e-04 0.002641350 4540  0.17222854   0.863
## conc0.12 B_2   1.792948e-03 0.002641350 4540  0.67880002   0.497
## conc0.4 B_2   -7.097261e-05 0.002641350 4540 -0.02686983   0.978
## conc1.2 B_2   -2.414977e-03 0.002641350 4540 -0.91429635   0.360
## conc4 B_2      4.571940e-03 0.002893454 4540  1.58009771   0.114
## conc12 B_2     3.670560e-03 0.002893454 4540  1.26857395   0.204
## conc40 B_2     3.970022e-03 0.002800092 4540  1.41781845   0.156

Significant differences are seen at 0.012 uM, 12 uM and 40 uM for B_0, and at 40 uM for B_1. Note that Loperamide was found active in Average Speed in Light and Average Acceleration in Light endpoints.

Evaluate t-tests for Dark.

##                       Value   Std.Error   DF     t-value p-value
## conc0.004 B_0 -0.6819917309 0.255309784  230 -2.67123225   0.008
## conc0.012 B_0 -0.3163190723 0.255309784  230 -1.23896181   0.216
## conc0.04 B_0  -0.3454960369 0.255309784  230 -1.35324245   0.177
## conc0.12 B_0  -0.3848596192 0.255309784  230 -1.50742213   0.133
## conc0.4 B_0   -0.2544640708 0.255309784  230 -0.99668750   0.319
## conc1.2 B_0   -0.4741939571 0.255309784  230 -1.85732779   0.064
## conc4 B_0     -0.4616918372 0.279677856  230 -1.65079869   0.100
## conc12 B_0    -0.3170755720 0.279677856  230 -1.13371711   0.258
## conc40 B_0     0.3552755277 0.270653603  230  1.31265767   0.190
## conc0.004 B_1 -0.0212202174 0.015956275 4540 -1.32989794   0.183
## conc0.012 B_1 -0.0294969973 0.015956275 4540 -1.84861423   0.064
## conc0.04 B_1  -0.0270402228 0.015956275 4540 -1.69464505   0.090
## conc0.12 B_1  -0.0299918156 0.015956275 4540 -1.87962512   0.060
## conc0.4 B_1    0.0054154895 0.015956275 4540  0.33939559   0.734
## conc1.2 B_1   -0.0168901299 0.015956275 4540 -1.05852586   0.289
## conc4 B_1     -0.0009266074 0.017479224 4540 -0.05301193   0.957
## conc12 B_1     0.0029855426 0.017479224 4540  0.17080522   0.864
## conc40 B_1     0.0346750870 0.016915229 4540  2.04993310   0.040
## conc0.004 B_2  0.0054204438 0.002289826 4540  2.36718574   0.017
## conc0.012 B_2  0.0030671315 0.002289826 4540  1.33946040   0.180
## conc0.04 B_2   0.0018925704 0.002289826 4540  0.82651271   0.408
## conc0.12 B_2   0.0001902576 0.002289826 4540  0.08308825   0.933
## conc0.4 B_2   -0.0002954939 0.002289826 4540 -0.12904644   0.897
## conc1.2 B_2    0.0014189362 0.002289826 4540  0.61966982   0.535
## conc4 B_2      0.0007913382 0.002508379 4540  0.31547792   0.752
## conc12 B_2     0.0010150525 0.002508379 4540  0.40466474   0.685
## conc40 B_2     0.0007595098 0.002427442 4540  0.31288482   0.754

Significant differences are seen at 0.004 uM for B_0, 40uM for B_1, and 0.004 uM for B_2. Note that Loperamide wasn’t found active in any endpoints in Dark.

Evaluate Phase Starting Points (Freeze and Startle)

In order to evaluate the starting points of each phase we need to find the estimate associated with each optimal model at t=0 (rather shifted t=-9.5) and extract a measure or error around this estimate in order to construct confidence intervals.

Start with Light model. Confidence intervals will be calculated about the fitted value at t=-9.5. SEs of these estimates will be calculated using the variance-covariance matrix of the fixed effects parameters. CIs will use t-tests with n-3 DF, 3 parameters per concentration, to assess accuracy of fitted values.

Calculate for treatment groups first.

fixedEf.L <- lmm2.L.REML$coefficients$fixed # Fixed effect parameter values
varFix.L <- lmm2.L.REML$varFix # VarCov matrix for fixed effects

y.hat.t_0 <- lapply(2:10, function(i) {
  
  t_0 <- -9.5 # Initial time point for centered data.
  
  # Extract parameter values to find model value at t = t_0
  B_0 <- fixedEf.L[i]
  B_1 <- fixedEf.L[11+i]
  B_2 <- fixedEf.L[20+i]
  
  # Extract Variance and Covariance of fixed effects to calculate SE of fitted value
  var.B_0 <- varFix.L[i,i]
  var.B_1 <- varFix.L[11+i,11+i]
  var.B_2 <- varFix.L[20+i,20+i]
  cov.B_0.B_1 <- varFix.L[i,11+i]
  cov.B_0.B_2 <- varFix.L[i,20+i]
  cov.B_1.B_2 <- varFix.L[11+i,20+i]
  
  # Find fitted value and SE
  mu.t_0 <- B_0 + B_1*t_0 + B_2*(t_0^2)
  se.t_0 <- var.B_0 + (t_0^2)*var.B_1 + (t_0^4)*var.B_2 + 2*(t_0*cov.B_0.B_1 + (t_0^2)*cov.B_0.B_2 + (t_0^3)*cov.B_1.B_2)
  
  # Calculate CI
  n = length( unique(Light.sqrt$fishID[Light.sqrt$conc==concs[i]]) ) # Number of observations
  alpha = 0.05
  CI.t_0 <- se.t_0*qt(1-(alpha/2), n-3)
  data.frame(mu.t_0 = mu.t_0,
             se.t_0 = se.t_0,
             CI.t_0 = CI.t_0)
})

y.hat.t_0.df <- cbind(conc=concs[-1], do.call('rbind', y.hat.t_0))

Calculate for control group.

  t_0 <- -9.5
  
  # Extract parameter values to find control model value at t = t_0
  B_0 <- fixedEf.L[1]
  B_1 <- fixedEf.L[11]
  B_2 <- fixedEf.L[13]
  
  # Extract Variance and Covariance of fixed effects to calculate SE of fitted value
  var.B_0 <- varFix.L[1,1]
  var.B_1 <- varFix.L[11,11]
  var.B_2 <- varFix.L[12,12]
  cov.B_0.B_1 <- varFix.L[1,11]
  cov.B_0.B_2 <- varFix.L[1,12]
  cov.B_1.B_2 <- varFix.L[11,12]
  
  # Find fitted value and SE
  mu.t_0.ctrl.L <- B_0 + B_1*t_0 + B_2*(t_0^2)
  se.t_0.ctrl.L <- var.B_0 + (t_0^2)*var.B_1 + (t_0^4)*var.B_2 + 2*(t_0*cov.B_0.B_1 + (t_0^2)*cov.B_0.B_2 + (t_0^3)*cov.B_1.B_2)
  
  # Calculate CI
  n = length( unique(Light.sqrt$fishID[Light.sqrt$conc==concs[1]]) ) # Number of observations
  alpha = 0.05
  CI.t_0.ctrl.L <- se.t_0.ctrl.L*qt(1-(alpha/2), n-3)
  
y.hat.t_0.ctrl.L <- c(mu.t_0=mu.t_0.ctrl.L, se.t_0=se.t_0.ctrl.L, CI.t_0=CI.t_0.ctrl.L)

Plot confidence intervals.

ggplot(y.hat.t_0.df, aes(y=mu.t_0,x=log10(conc))) +
  geom_point() +
  geom_errorbar(aes(ymin=(mu.t_0-CI.t_0), ymax=(mu.t_0+CI.t_0))) +
  geom_hline(yintercept=0) +
  labs(title = "Light LME Model: Fitted Value at Initial Time interval",
       subtitle = "Intersection of CI with y=0 indicates insignificant t-test p-value",
       y = "Fitted Value - (Control Fitted Value)")

Repeat for Dark model. Code for calculating confidence intervals won’t be shown, but method is identical to method above.